Se pretende estimar \(\mathbb{E}_{f}[g(X_{1}, X_{2})]\), donde \[ f(x_{1}, x_{2}) \propto \begin{cases} \mathrm{e}^{-4 (x_{2} - x_{1}^{2})^{2} + (x_{2} - 1)^{2}} &\text{si } x_{2} \leq 2, \\ 0 &\text{en otro caso.} \end{cases} \qquad g(x_{1}, x_2) = \sqrt{x_{1}^{2} + x_{2}^{2}} \]

library(plotly)

f_u <- function(x1, x2) {# fdd no normalizada
  exp(-4 * (x2 - x1 ^ 2) ^ 2 + (x2 - 1) ^ 2) * 
    (x2 <= 2) # multiplico por la condición para que fuera una función vectorizada
}

# malla de puntos donde voy a representar la función
malla <- data.frame(x1 = (-250:250) / 100,
                    x2 = (-250:250) / 100)
plot_ly(malla, x = ~x2, y = ~x1) %>% 
    add_surface(z = outer(malla$x1, malla$x2, f_u) # outer: aplique la función en todas las combinaciones de la malla
                )

Tiene una zona donde acumula mayor probabilidad y dos picos donde hay menor probabilidad.

Vamos a realizar la estimación requerida a partir de un paseo aleatorio de Metropolis, haciendo uso del paquete mcmc. Para ello, en primer lugar debemos definir el logaritmo de la densidad objetivo como una función que toma como argumento un vector con los dos valores \(x_{1}\) y \(x_{2}\).

log_f_u <- function(x) { # función de un único argumento 
                         # (vector del estado, el vector (x1,x2))
  x1 <- x[1]
  x2 <- x[2]
  if (x2 <= 2) {
    -4 * (x2 - x1 ^ 2) ^ 2 + (x2 - 1) ^ 2
  } else {
    -Inf
  }
}

A continuación, iniciamos el paseo aleatorio de Metropolis a partir del estado \(x_{1} = x_{2} = 0\), que según el gráfico se encuentra cerca de donde la densidad toma su máximo (donde acumula mayor probabilidad). Podemos comprobar que la función metrop devuelve la realización de la cadena de Markov en forma de matriz con dos columnas, ya que cada estado es una tupla de dos números reales.

Hemos tomado el inicial (0,0) porque en el gráfico vemos, mas o menos, que en el (0,0), nuestra curva acumula la mayor probabilidad. Hay que escoger un punto en la zona de probabilidad, aproximado.

paseo_Metropolis <- mcmc::metrop(log_f_u, initial = c(0, 0), nbatch = 1e4)
head(paseo_Metropolis$batch)
##           [,1]       [,2]
## [1,] 0.0000000  0.0000000
## [2,] 0.0000000  0.0000000
## [3,] 0.0000000  0.0000000
## [4,] 0.7588121 -0.2884623
## [5,] 0.5888510 -0.2832625
## [6,] 0.5888510 -0.2832625

Rechazo el primero tres veces, luego acepto luego acepto y luego rechazo (en x1). Para x2, rechazo tres veces y luego acepto y vuelvo a rechazar dos veces.

La regla empírica establece que para acelerar la convergencia de la cadena de Markov a la densidad objetivo, el porcentaje de estados aceptados debería ser del 25 %.

paseo_Metropolis$accept # pordríamos admitir, está cerca del 25%
## [1] 0.2293
# Para acercarme al 25% de estados aceptados, disminuyo la varianza
paseo_Metropolis <- mcmc::metrop(paseo_Metropolis, scale = .9) 
# scale: controls the proposal step size
paseo_Metropolis$accept # me he pasado, voy a aumentar la varianza un poco para disminutir el porcentaje de estados aceptados
## [1] 0.2662
paseo_Metropolis <- mcmc::metrop(paseo_Metropolis, scale = .95)
paseo_Metropolis$accept
## [1] 0.2416
paseo_Metropolis <- mcmc::metrop(paseo_Metropolis, scale = .93)
paseo_Metropolis$accept
## [1] 0.2471

Como diagnóstico de convergencia, los siguientes gráficos parecen indicar que la cadena de Markov ha convergido a la densidad deseada.

library(coda)
traceplot(mcmc(paseo_Metropolis$batch))

# Me da dos gráficos, porque estoy en R2, uno para cada dimensión, x1 y x2 respectivamente
library(ggplot2)



ggplot(expand.grid(# expand grid. Hace todas las combinaciones, es una malla
  x1 = (-250:250) / 100, # lo que hago es: seq(-2.5 , 2.5 ,by=0.01)
  x2 = (-250:250) / 100), # lo mismo
       aes(x = x1, y = x2, z = f_u(x1, x2))) +
  geom_contour() + # gráfico de líneas de nivel
  geom_point(data = as.data.frame(paseo_Metropolis$batch), # puntos a partir de los estados
             mapping = aes(x = V1, y = V2),
             size = .5,
             inherit.aes = FALSE)

Con respecto a las dos primeras gráficas:

La estabilización de los valores en un rango en principio mostraría convergencia, pero el gráfico no es absoluto, en realidad habría que probarlo matemáticamente.

La cadena de markow ha convergido a la zona de probabilidad.

En el último gráfico veremos las curvas de nivel, que juntarán puntos que tienen el mismo valor de la función de densidad. Si están muy juntas la fdd tiene mucha pendiente y si se separan son zonas donde la fdd tiene menos pendiente. Las curvas de nivel me dan los tres picos de la fu, los puntos se situan en mayor medida donde están esos tres picos, mayoritariamente een el centro.

Hay puntos fuera poque la probabilidad ahí no es nula, existe probabilidad de que se generen esos estados.

Estamos ya en condiciones de obtener la estimación requerida, así como un intervalo de confianza para la misma. El gráfico de autocorrelación parece indicar que, al aplicar el método de las medias por lotes, la longitud de estos debería ser de al menos 50 estados.

autocorr.plot(mcmc(paseo_Metropolis$batch), lag.max = 100)

gráfico de autocorrelación para x1 y x2:

Por este motivo, se ve que a partir de 150, las dos coordenadas son independientes, por tanto, debo coger una longitud de lotes de 150

g <- function(x) {
  x1 <- x[1]
  x2 <- x[2]
  sqrt(x1 ^ 2 + x2 ^ 2)
}
paseo_Metropolis <- mcmc::metrop(paseo_Metropolis, outfun = g)
autocorr.plot(mcmc(paseo_Metropolis$batch), lag.max = 100)

# aptes día 7/12
library(mcmc)

# puedo decirle a metrop el número de lotes y el número de estados en cada lote, ahora necesito decirle que antes de hacer la media aplique g a cada estado(de 1 a 10000), ya que yo quiero la media de g(estados)
paseo_Metropolis <- mcmc::metrop(paseo_Metropolis, nbatch = 1e4, blen = 1e2, 
                                 outfun = g) # función de salida (recibido en un vector)
autocorr.plot(mcmc(paseo_Metropolis$batch), lag.max = 100) # gráf de autocorrelac de x1 y x2

# autocorrelación cuando YA HE APLICADO G

Cuando ya he aplicado g, solo obtengo una gráfica, ya que por definición de la función \(g(x_1,x_2)\), obtengo un único número para cada estado. Alcanzamos la independencia en 40 y pico, pero elijo longitud 100 para asegurar la independencia.

valores <- as.vector(paseo_Metropolis$batch)
estimacion <- mean(valores)
error_estandar <- sqrt(var(valores) / paseo_Metropolis$nbatch)
probabilidad_cobertura <- 0.95
alfa <- 1 - probabilidad_cobertura
percentil <- qnorm(1 - alfa / 2)
intervalo_confianza <- estimacion + c(-1, 1) * error_estandar * percentil

Se obtiene entonces una estimación de 0.864163, con (0.8582931, 0.870033) un intervalo de confianza con probabilidad de cobertura 0.95.